Pseudogap and preformed pairs in the imbalanced Fermi gas in two dimensions 



<Z3 



S.N. Klimin, 1 ^ J. Tempere, 1 ' 2 and J.T. Devreese 1 '^ 

1 Theorie van Kwantumsystemen en Complexe Systemen, 
Universiteit Antwerpen, Universiteitsplem 1, B-2610 Antwerpen, Belgium 
2 Lyman Laboratory of Physics, Harvard University, Cambridge, MA 02138, USA 

(Dated: July 12, 2012) 

The physics of the pseudogap state is intimately linked with the pairing mechanism that gives 
1 rise to superfluidity in quantum gases and to superconductivity in high-T c cuprates, and therefore, 

both in quantum gases and superconductors, the pseudogap state and preformed pairs have been 
under intensive experimental scrutiny. Here, we develop a path integral treatment that provides a 
£Nj 1 divergence-free description of the paired state in two-dimensional Fermi gases. Within this formal- 

ism, we derive the pseudogap temperature and the pair fluctuation spectral function, and compare 
these results with the recent experimental measument of the pairing in the two-dimensional Fermi 
gas. The removal of the infrared divergence in the number equations is shown both numerically and 
analytically, through a study of the long-wavelength and low-energy limit of the pair fluctuation 
density. Besides the pseudogap temperature, also the pair formation temperature and the criti- 
cal temperature for superfluidity are derived. The latter corresponds to the Berezinski-Kosterlitz- 
Thouless (BKT) temperature. The pseudogap temperature, which coincides with the pair formation 
temperature in mean field, is found to be suppressed with respect to the pair formation temperature 
by fluctuations. This suppression is strongest for large binding energies of the pairs. Finally, we 
investigate how the pair formation temperature, the pseudogap temperature and the BKT temper- 
ature behave as a function of both binding energy and imbalance between the pairing partners in 
the Fermi gas. This allows to set up phase diagrams for the two-dimensional Fermi gas, in which 
the superfluid phase, the phase-fluctuating quasicondensate, and the normal state can be identified. 
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^3 ■ I. INTRODUCTION 

a: 

Ultracold atomic gases are increasingly used as quantum simulators to probe many-body physics Recent 
efforts have focused in particular on understanding superfluidity and Cooper pairing in interacting Fermi systems, 
including i.a. the effects of varying the interaction strength, introducing population imbalance, and reducing the 
dimensionality. When these fermionic superfluids are described in the path integral formalism, the thermodynamic 
potential ft(T, fif, /ijj of the interacting Fermi gas (as a function of temperature T and chemical potentials m,(i\, of 
spin-up and spin-down components) is rewritten as a functional integral over a bosonic field A XiT that represents the 
field of the pairs. This field is introduced through the Hubbard-Stratonovich transformation that allows the exact 
f**^ " elimination of the fermionic degrees of freedom and results in an action functional for the bosonic field 0. At first sight, 

j. | , one has only succeeded in rewriting an unsolvable functional integral over fermionic fields by an equally unsolvable 

functional integral over the bosonic fields. However, the bosonic field lends itself to an obvious simplification when one 
j intuits that a (uniform) Bose-Einstein condensation of pairs is present. In that case, one can surmise that A x r « A, 
7-H ■ i.e. all pairs are in the zero-momentum state so the field is a constant in real space. The functional integral can 
then replaced by its saddle-point value, substituting A X)T ~ A and dropping the integrations. The optimal value of 
A is found by extremizing the action or, equivalently, by minimizing the thermodynamic potential fl sp (T, fj^, fj,±; A) 
with respect to the saddle point. When A = 0, the (non-interacting) normal Fermi gas is obtained, when A^O, the 
J-j . saddle-point approximation bears out a Bose-Einstein condensate of pairs. 

In two-dimensional superfluid systems, the interpretation of the bosonic field is more subtle. To be precise, the 
condition for Bose-Einstein condensation (BEC) has been defined by Penrose, Onsager [5| and Yang [6[ as the presence 
of off-diagonal long range order, i.e. linix-x/^oo (A X)T A X / )T ) =^0. At the level of the saddle point, A XjT ss A, it is clear 
that a nonzero A implies BEC. However, as Mermin,Wagner Q and Hohenberg 0] pointed out, in the two-dimensional 
system fluctuations play a crucial role: they will prohibit off-diagonal long range order in uniform systems. These 
fluctuations around the saddle point are commonly taken into account through the Bogoliubov shift A X T = A + </> XiT , 



o 
o 



> 
o 



a On leave of absence from:Department of Theoretical Physics, State University of Moldova, str. A. Mateevici 60, MD-2009 Kishinev, 
Republic of Moldova. 

b Also at:Technische Universiteit Eindhoven, P. B. 513, 5600 MB Eindhoven, The Netherlands. 



2 



whereafter x r is treated as a small fluctuation so that only terms up to second order in XjT are retained in the action 
functional. Then, for a given A, the thermodynamic potential f2(T, A) is expressed as a functional integral 

over the fluctuation fields, with a quadratic action. The quadratic functional integral can be performed, and we obtain 
a fluctuation correction Qfi(T, /i^, fj,y, A) = Q(T, /i-j-, fi^; A) — Q sp (T, /i-j-, fj,^; A) to the thermodynamic potential. The 
fluctuation fields need not be written down as complex fields x , r resulting from the Bogoliubov shift: equivalent 
results are obtained by introducing (real) amplitude and phase fluctuation fields through A XjT ss A (1 4- S X . T ) e l8x - T . 

The culprit suppressing off-diagonal long range order in 2D is precisely the phase fluctuation field e t6 ^- T . Indeed, 
Mermin and Wagner show that (A X)T A X / )T ) « /A 2 e* ex T ) — > due to the long-wavelength behavior of X . T , the relative 
phase. According to the Penrose-Onsager-Yang criterion this means that Bose-Einstein condensation does not occur. 
However, we can identify other interesting phases from a study of the bosonic pair field A x T . Firstly, {V 0X ' T ) — > 
does not imply that A = 0, as noted by Kagan in his study of quasicondensation 0- We can identify A ^ with 
the presence of pairing, and search for a transition temperature T* for pair formation separating the A = phase 
from the A ^ low temperature phase. Second, although BEC is suppressed, superfluidity can still be present in 
the two dimensional system below the Berezinski-Kosterlitz-Thouless [H, Q (BKT) temperature Tbkt- The order 
parameter for superfluidity is p s , the superfluid density, defined as the phase stiffness and calculated as the prefactor 
of the (V# XjT ) 2 term in the Lagrangian for the phase field, as explained in more detail below. Kosterlitz and Thouless 
Q described a mechanism whereby phase stiffness can be lost, namely through the appearance and unbinding of 
vortex-antivortex pairs that start to proliferate at Tbkt and scramble the phase field. This mechanism was observed 
experimentally in an 2D atomic Bose gas by Dalibard and co-workers (Io| . 

When the bosons under consideration are composite particles, such as Cooper pairs, a third relevant temperature can 
be identified, related to the density of states of excitations, or equivalently the spectral function for the fluctuations. 
As we will show below, the fluctuation terms in the density can be expressed through a spectral function g(q, uj) 
describing the contribution of fluctuations with a given wave number q and momentum ui. In the Nozieres and 
Schmitt-Rink (NSR) formalism for the two dimensional system [l2|], the integral over the spectral function is 
divergent, invalidating the number equations = n a , a =f, i- We show that this divergency is absent when we 

apply the formalism of Hu, Liu and Drummond [13l4l5| , which these authors dubbed the Gaussian Pair Fluctuation 



(GPF) approach, to the two dimensional case. As we show below (section |TTC|, the GPF approach allows to set up and 
simultaneously solve the gap and number equations also in the two-dimensional case. This allows us to derive results 
for A,p s that take into account fluctuations (both phase fluctuations and amplitude fluctuations). The resulting 
fluctuation spectra can then be used to obtain the finite-temperature thermodynamics of the two-dimensional Fermi 
superfluid, following the approach of Salasnich for the three-dimensional case [16| . 

In studying the fluctuation spectra, we find that for temperatures above a critical temperature T pi the fluctuation 
spectral function <?(q, w) becomes negative at long wavelengths (q < q c ). This happens at a temperature above Tbkt 
and (obviously) below T*. We interpret this temperature T p as the pairing temperature at which the pseudogap is 
open, inspired by the recent experiments (T3, EH which investigated pairing in ultracold 2D atomic Fermi gases. For 
the ultracold atomic gases in 3D, the pseudogap state above the critical temperature is a subject of an intense study, 
both experimental and theoretical |19l424ll . and the similarity with the pseudogap physics in superconductors has not 
gone unnoticed [25] . In the experiment [l 7j j , the spectral function for excitations of the Fermi gas is determined through 
momentum- resolved photoemission spectroscopy. In Ref. jl8|, the momentum-integrated photoemission spectra are 
measured, and the evolution of fermion pairing was followed from three to two dimensions by varying the strength 
of the confining optical lattice. Both experiments reveal a non-zero pairing gap. While the experiments in [17| were 
interpreted to reveal the pseudogap, i.e. pairing in a non-superfluid state, superfluidity itself has not been observed 
in 2D yet. Therefore, the existence of the pseudo-gap regime is not experimentally settled until superfluidity itself is 
observed at temperatures lower than the temperature for pair formation. 

We compare the fluctuation spectral functions derived from our microscopic (GPF-based) theory to the measured 
spectral functions for excitations, and we also compare the measured pseudogap temperatures with the calculated 
T p as a function of the interaction strength. We find that fluctuations indeed greatly lower the temperature range 
of existence of the pseudogap phase, especially in the strong-coupling regime. Consequently, in order to obtain a 
complete phase diagram for the Fermi gas in 2D, we must consider each phase taking in account fluctuations. To 
the best of our knowledge, this problem has hitherto not yet been satisfactorily solved for the case of fermions in 
2D because of the aforesaid divergence of the density due to fluctuations at finite temperatures. Here, as mentioned 
above, we tackle the problem by correcting the NSR approach using the Gaussian Pair Fluctuation theory (GPF) 
proposed by Hu, Liu and Drummond |13l4l5| for the three dimensional Fermi gas. Moreover, we extend the results 
to the case of imbalance. 




The paper is organized as follows. In Sec. HU we present the divergence- free method for the self-consistent calculation 
of thermodynamic parameters of interacting imbalanced fermions in 2D taking into account both amplitude and phase 
fluctuations. In Sec. IIIIl density distribution functions for an imbalanced 2D Fermi gas are investigated. In Sec. IIV1 
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we discuss finite-temperature phase diagrams for the imbalanced Fermi gas in 2D. In Sec. |Vj the theory is applied to 
the interpretation of the experiment on pairing of cold atoms in 2D. The discussion is followed by conclusions, Sec. 

m 

Before going forward in the next section with presenting the functional integral approach in the GPF framework, 
it is useful to note that the GPF approach that we follow here is not the only way to avoid the divergence problem 
that occurs in the NSR description for the Fermi gas in two dimensions. The NSR scheme and its modifications are 
related to the T-matrix perturbation approach, in which the effective interaction between pairs is taken into account 
diagrammatically. In this context, the divergence problem for a Fermi gas in two dimensions can be remedied by 
taking into account higher orders of the T-matrix expansion - via an effective interaction between pair fluctuations 
[29| . This interaction stabilizes the superfluid phase of the 2D fermion system at very low temperatures. However, 
in 2D the T-matrix method does not predict the universal jump in the superfluid density (30j related to the BKT 
phase transition. A correct description of the superfluid density becomes possible by explicitly focusing on the phase 
fluctuations, as in the approach of Refs. (26U281 . l3lj|. Within that approach, bosonic pair field is gauge transformed 
A x , r e I03C ' T and a subsequent gradient expansion of the fluctuation action is performed for phase fluctuations assuming 
that phase gradients are small. This leads to a quadratic effective action functional of the phase field 9, which, as 
distinct from the scheme of Ref. [Ill], has no divergence for A^O. As far as gradients of the fields are assumed to be 
small, the resulting effective action is treated as a hydrodynamic action (see, e.g., Ref. [Hj]). The gradient expansion 
does not contain the a priori assumption that fluctuations themselves are small. In this connection, the method was 
categorized in Ref. [26] as non-perturbative. In Ref. (33[, the present authors applied the method of Refs. [26T - I281 l3l| 
to derive the effective hydrodynamic action for a Fermi gas with a population imbalance. A non-perturbative approach 



was also the key to develop a description free of infrared and ultraviolet divergences for the 2D Bose gas |34[, that 
successfully describes the crossover between the mean-field regime and the critical fluctuation range corresponding 
the BKT transition HI. 



II. THERMODYNAMIC FUNCTIONS OF THE FERMI GAS IN 2D 



A. Gap equation 

We consider a gas of interacting fermions in 2D, with a contact interaction and with s-wave pairing. In the ultracold 
regime where only s-wave interactions matter, these interactions only take place between "spin-up" and "spin-down" 
fermions (in practice, these are usually two different hyperfine states of an atomic species). The thermodynamic 
functions of the Fermi gas are completely determined by the partition function. Here we will focus on the thermody- 
namic potential fi per unit area. The treatment is performed within the path-integral formalism following Ref. [33j j . 
building on the original path-integral treatment in Ref. @ for the case of a balanced three-dimensional Fermi gas. 
The partition function is represented as the path integral over Grassmann variables ipa- (x, r) , ipg. (x, r), 

Z = e-WT'KQ = ( £^,x,T^a,x,T exp (-S) . (1) 



The action functional of interacting fermions is given by the integral 



S = J dr J d 2 x ^2 ^o-,x,t (J^ - V x - fi a ^j ^<r,x,r 



+ g I dr I d 2 x V't, x ,rV'4.,x,rV'4.,x,T^t,x,T, (2) 



' 2, 







where g is the interaction strength and (3 — l/(fcgT) is the inverse thermal energy. We choose a system of units where 

1/2 

H = 1, 2m = 1, and the Fermi wave vector kp = (27m) = 1 with n the total density. Here, we consider also the 
case when imbalance is present, i.e. the number of spin- up and spin-down atoms are unequal: ^ n^. This in turn 
implies that the chemical potentials m and m should be fixed separately. Rather than contemplating the separate 
components, we will work with the total density n = n^ + nj, and the density difference Sn = n^ — n^. Correspondingly, 
we will use the average chemical potential /1 = (/i-j- + //j,)/2, and the chemical potential difference £ = (/i-j- — /zj,)/2. 
Note that the total density is equal to l/(27r) in our units, so this means that we need to solve the number equation 
to fix [i (in the non-interacting case, fi — 1 in our units). Only with respect to the imbalance, we have a choice of 
studying the free energy (and making phase diagrams as a function of cm) or the thermodynamic potential (and 
making phase diagrams as a function of £). The thermodynamic potential is linked to the free energy by the usual 
Legendre transform and, as mentioned, in our formalism this corresponds to imposing the number equations. 
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The strength g of the contact interaction is renormalized as in Refs. 3G, 37] using the binding energy Ei, for a 
two-particle bound state, which always exists in 2D |43l |44|| : 

i = jLf ]B * +< ;|-/jgfc . 1 _ ( 3) 

9 &tt V E ) J (2nf2k 2 -E + i5 y ' 

with 6 a positive infinitesimal number. The BCS regime corresponds to E^/Ep <C 1, whereas the BEC regime 
corresponds to the opposite ratio Eb/Ep ^> 1. Similarly to Ref. Q, we introduce the pair field A x r and perform 
the Hubbard-Stratonovich transformation, which results in a fermion-boson action quadratic in fermion variables. 
After integrating out the fermion variables, the following effective bosonic action is obtained as a functional of the 
Hubbard-Stratonovich pair field A XiT : 



S eff = - tr [In - dr J d 



2 A X , T A : 



x ^" x ' T , (4) 



where G 1 is the inverse of the Nambu propagator 



= CT o ( ^ -C) -ois (V 2 -<tiA x , t . (5) 



Here, aj are the Pauli matrices. As far as the effective action S e ff is not a quadratic functional of the Hubbard- 
Stratonovich pair field, the resulting functional integral over the pair field 

Z oc J VA X:T VA x , r exp(-S ef f) (6) 

cannot be calculated analytically exactly. As in the analogous problem in 3D 0, [3|1 HH , and as explained in the 
introduction, we consider approximations provided by an expansion of the effective action S e jf over fluctuations 
of the pair field A x T about its saddle-point value A. The phase diagrams of a 2D Fermi gas in the saddle-point 
approximation have been investigated in Refs. [37ll40l|. The effective saddle-point action provides the thermodynamic 
potential per unit area: 

Q sp (T, M ,C;A) = - 



d 2 k 

(2^f 



In (2 cosh /3E k + 2 cosh /3Q 



A 



2 



(7 » 



Here, ^ = k 2 — /i is the fermion energy, and E^ = y/& + A 2 is the Bogoliubov excitation energy. The gap parameter 
A is determined from the gap equation generalized to the imbalance case - the minimum condition for the saddle-point 
thermodynamic potential as a function of the gap parameter A at fixed temperature and chemical potentials: 

dn sp Q3,M;A) 

dA = °' (8) 

For high temperatures (T > T*) or at high levels of imbalance (£ > Cc), thermodynamic potential will have its 
minimum at A = 0, the unpaired normal state. Following the experimental observation of superfluidity in imbalanced 
Fermi gases in 3 D I 4llj , the phase diagram of the imbalanced Fermi gas has attracted a lot of attention (for a recent 
review, see Ref. [42j). To find the phase diagrams in 2D, the above gap equation has to be solved in conjunction with 
the number equations discussed in the remainder of this section. 



B. Gaussian fluctuations 

The next-order approximation brings into account fluctuations about the saddle point: 

A XiT = A + XjT , 
A X , T = A + X , T . 

We apply the Fourier expansion of the fluctuation coordinates: 



1 oo 

4 - = 7^E E e^ r -^> k K), (9) 



q n— — oo 

oo 



- = 7^E E e - iq ' r+ ^ q K) (io) 



v ' q ri— — oo 
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where L is the linear size of the 2D system, and u> n = 2irn/ '(3 (with n = 0, ±1,±2, . . .) are the bosonic Matsubara 
frequencies. The quadratic fluctuation contribution to the effective bosonic action is the functional of complex 
fluctuation coordinates similar to that derived in Ref. 391: 



) ¥>q (W„) 



q n— — oo 



+<^ q (u; n )^_ q (u;_„)]}, (11) 

where Afj^ (q,iui n ) are the matrix elements of the inverse pair fluctuation propagator. They are determined by the 
expressions (cf. Ref. [Hj]): 

1 f d 2 k X(E k ) 



_ _ , . . i r d J ] 

M hl (q,iu] n ) = — + — 
9 J (2tt 



(iUJ n - g k + gk+q) (-Ek + gk) 

^iu;„ - £k + £k+q) («W n - -Ek - E k+q ) 

(j^n + E k + gk+q) (ffk - 60 

(iw n + £k - -Bk+q) (iw n + Ek+q + -E k ) 

and 



(12) 



Mi, a (q,»w„) = —A / — 

J (27T 



<Pk X(E k ) 
(2tt) 2 2£ k 

1 



(lW„ — i?k + -Ek+q) (i^n — Ek — i?k+q) 
1 



(iw„ + E k - E k+Cl ) (iuj n + E k + £k+q) 
Here, the following function has been introduced, 



(13) 



X (F\= smh(/3E k ) 

1 k; cosh(/3£ k ) + cosh^C) ' 

The integration over fluctuation coordinates gives us the fluctuation contribution f2^(T, /i, £; A) to the total grand- 
canonical thermodynamic potential f2 per unit area: 

fi i7 (T, m>C;A) = — / —^2 E ln F (<*< ^ ( 15 ) 

zp J (2tt) n= _ 0o 



with 



T (q, ift n ) = Mi,! (q, iw„) Mi,i (q, -iw n ) 

- Afi )2 (q, «w„) Mi )2 (q, -«w„) . (16) 



C. Number equations and the GPF approach 



The fermion density and the density difference fix the chemical potentials [i and £ through the derivatives of the 
total thermodynamic potential per unit area: 



an 



(17) 
(18) 
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(remember that in our units n = 1/2%). We can write out these equations by splitting the total thermodynamic 
potential in saddle point and fluctuation contributions. 

(19) 

Sn : — ^ (20) 



d£2 sp 








dp 


T,C 


dp 








dVL n 




dc 


T,n 


d( 





n f i, 



We will denote the first and second terms in the right hand side (RHS) of expression (fl9l) for n as n sp and 
respectively. Similarly, the terms in the RHS of expression (|20|) will be denoted by Sn sp and Srifi. Note that the 
thermodynamic potentials obtained from expressions (0 and (|15|) are expressed as a function not only of T, p, £, but 
also of A. This gap A is not an independent thermodynamic variable, and when considering f2(T, p, £, A) explicitly as 
a function of also A, the implicit dependence of A on the chemical potentials must be taken into account in (fT9|) . (|20l) : 



8n 













dA 




dfi 


TX-A 


dp 


TX-A 


dA 




TX 






dn fl 




dClfi 


dA 




dc 


T,n,A 


d( 


T,nA 


dA 


tx,» ~^ 





(21) 



Note that the gap equation dtt sp /dA — at fixed T,p,( implies that the implicit dependence of A on the chemical 
potentials will only affect the fluctuation part of the thermodynamic potential in the above equations. Different 
theories of the BEC-BCS crossover, in any dimension, can be categorized by their choice of number and gap equations. 
The simplest mean field approach only keeps the terms with Q sp . The Nozieres and Schmitt-Rink approach also 
includes the second terms in the RHS of expressions (|21[) . Finally, the Gaussian pair fluctuation approach includes 
also the last term in the RHS of expressions (|21l) . Note that in the literature, there is no common opinion on which 
approach is best. For example, on the one hand, Randeria et al. 0, Hfjj], Hu et al. [l3j], Keeling et al. [13] state that 
the derivative over p must be performed taking into account a variation of the gap determined by the gap equation. 
On the other hand, Ohashi et al. [48l - [50| , and Strinati et al. [Hi HH use the other definition, considering A in the 
number equations as an independent variable and, therefore, applying the gap equation after taking the derivatives 
dfl/dfi. In the papers |48l - [50l |. it is stated that the last terms in (|2"Tj) are the higher-order corrections with respect 
to Gaussian quadratic fluctuations. Keeling et al. (47j correctly argue that both terms in those derivatives are of 
one and the same order and emphasize that the existence of the second term is crucial in two dimensions. Below, 
we demonstrate the key significance of taking into account of the last terms ([2~Tj) for the convergence of fluctuation 
contributions to the fermion density in 2D. 

As stated in the introduction, in order to treat the fluctuations, there is an alternative to the Bogoliubov shift 
A XjT = A + </>x,r, namely the parametrization in amplitude and phase fluctuations A X . T ~ A (1 + 8 X . T ) e l6 *- T . When, 
after this parametrization, the effective action is expanded with respect to <5 Xi7 - and 9 X _ T (rather than XjT and </> x ,t), 
this leads to another quadratic fluctuation action S'^ which differs from expression (jlll) for 5/; only by terms that 
vanish when applying the gap equation. Correspondingly, the thermodynamic potentials O/; and £1'^ provided by 
those two actions lead to one and the same contribution to the fermion density. Furthermore, keeping in S'^ only 
the leading order long-wavelength and low-energy terms leads to the same effective "hydrodynamic" action as in Ref. 
[33} . In the particular case of a balanced gas, the effective action of Ref. [33| turns to the result of Refs. [2(1 l3lj . 
This means that the effective action described as the result of the non-perturbative approach in Ref. (26[ can be 
equivalently re-derived within the perturbative NSR-like scheme. Moreover, the present treatment can be considered 
as an extension of the approach of Refs. [2(| HH HH beyond the long-wavelength and low-energy approximation 
(and to imbalanced 2D gases). The hydrodynamic action is particularly useful in extracting the superfluid density 
p s , by identifying it with the prefactor of the (V0 Xir ) 2 /2 term in the expression for S'^. This identification yields 
straightforwardly^ : 



p a (T,n,C,A) 



1 

4-7T 



-II.- I< I 1 - ^rX(Ek) - k 2 X'(E k 



(22) 



with X(Ek) given by expression (fT4|) and X'(Ek) its first derivative, evaluated in Ek- Once A,//, £ are obtained for 
a given temperature (and interaction strength) by solving the gap and number equations, they can be substituted in 
this expression to determine whether the system is in the superfluid phase (p s 7^ 0) or the normal phase (p s = 0). 
As discussed in the results section, we also use this expression to find the temperature Tbkt of the phase transition 
between those two states. Already we note that A = leads to p s — 0, so that Tbkt < T* and the superfluid state 
requires pair formation, as it should. 
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D. Pair fluctuation spectral functions 



From expression (J7J for the saddle-point thermodynamic potential, we derive the following expressions for the 
saddle-point densities: 



an. 



dn sp := — 



d 2 k 

d 2 k 



1 



sinh(/3^ k ) 



£ k cosh (/?C) + cosh (PEy,) 
sinh (PC) 



(2tt) 2 cosh (/3C) + cosh (££ k ) ' 



(23) 
(24) 



Similarly, the fluctuation contributions to the fermion densities are determined using (|f 51) and using the GPF approach. 
The results can be written as a sum over wavelengths and frequencies of pair fluctuation structure factors: 



on 



n f i := - 



Srifi := 



n 



dn fl 



T,ii 



/lRJl J(w) 

r d 2 q 1 ^ 



(25) 
(26) 



The pair fluctuation structure factors J and K are given by 



J (q,iu; n ) 



if (q, to; n ) 



1 



5M 1; i (q, iw r , 



T(q, iu„) 
9M ij2 (q, iw. 



1 



-Mi, 2 (q, -iw n ) 



Mi,i (q, -h) 



T(q, iw n ) 



9Mi,i (q, iw n ) 



Mi,i (q, -««; n ) 



9Mi, 2 (q, MJn) 



Mi, 2 (q, -iw„) 



(27) 



(28) 



We transform the Matsubara summations in (j2"5)l and (|26l) to the contour integrals in the complex plane as follows: 

Im J (q, cj + i£) 



1 1 / 

- V J (q, iw„) = - / 



1 



-dio, 5 -> +0. 



(29) 



This allows to express the resulting fluctuation contributions to the fermion density through the distribution functions 
for pair excitations: 



rifi 
Srifi 



1 f°° 
= T~2 / 9n (q) qdq, 

gsn (?) grfg- 



1 

2^2 



(30) 
(31) 



The fluctuation distribution functions g n (q) and gg n (q) are the integrals over the frequency with the pair fluctuation 
structure factors: 



9n (?) = - 
9Sn (q) = - 



Im J (q, u> + iS) 

eP u - 1 
Imif (g,u; + i<5) 



1 



(32) 
(33) 



The functions g n (q) and gs n (q) are proportional to the densities of states for the pair fluctuations. The behavior 
of these functions is crucial for understanding of the pseudogap properties and of different phase transitions in the 
imbalanced 2D Fermi gas. In the NSR approach, the fluctuation distribution functions have a divergency, and as a 
consequence no value of the chemical potential /i can be found so that the number equation n sp + rifi = n = 1/ (2n) 
is satisfied. In the GPF approach, the divergency is overcome and the number equation can be satisfied. In order to 
demonstrate this, we focus in the next section on the long wavelength limit where exact analytic expressions for the 
distribution functions are obtained. 
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III. DISTRIBUTION FUNCTIONS IN THE LONG- WAVELENGTH LIMIT 

A. Long-wavelength expansion 

In order to investigate the problem of the long-wavelength convergence for the fluctuation contributions to the 
density, it is necessary to derive analytically the spectrum of low-lying and long-wavelength pair excitations. For this 
purpose, we expand the matrix elements of the inverse pair fluctuation propagator Mjk (q, z) in powers of (q, z) up 
to the second-order terms in power of z and q, 

M M (q,z)« A + Bq 2 + Cz + Fz 2 , (34) 
Mi, 2 (q,z)w D + Eq 2 +Hz 2 . (35) 

The derivation of the coefficients is rather tedious. In this connection, here only the final results are represented. The 
coefficients in the Mxi expansion are given by 

A -±_ r kdk ( 1 E ' + ^X(E„) g*Wl (36) 

B = j_ r kdk ^{Ei + iE 2 e k -mp-E^ + nEi x 

167T J El 1 

A 2 [°° h ,J §k (g fc 2 -3fc 2 6) 
+ ^J o kdk^ ^ X(E k ) 

+ 6(3^-|)-^ ^ (£fc) _^ (3)( ^ i (37) 
C-ijT^W-^jT**^ (38) 

F =-wJ kdk ^ET x ^- (39) 

and those in the M%2 expansion are given by 



a 2 r ,^ iok2 e k -3E 2 ^ k +k 



2\ 



E=— kdk ^ -fAii >- X (£ k ) 

16tt ,/o 



57T 

-^ (3) (i? fe )), (41) 



2E e k 



32ttJ q El 

In these expressions X',X" and X^ are the first, second and third derivatives of the function X given by expression 
(|14p . with respect to its argument. In the literature, an analogous expansion was performed for 3D at finite temper- 
atures in Refs. p8l - l50l | in the strong-coupling limit, and in Ref. [45| at low temperatures. The present expansion is 
all-coupling and all-temperature, because no restriction is imposed on the thermodynamic parameters. 
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B. Structure factor 

Let us substitute the expansions (|34[) . (|35|) to the structure factor J (q, z) in order to obtain its long- wavelength 
and low-energy form Ji w (q, z). The spectrum of pair bosonic excitations is determined by the poles of Ji w (q, z). 
In the long-wavelength and low-energy range, these roots are z = ±u> q with the pair excitation frequency u q which 
satisfies the Goldstone theorem, 



cu q = q\Jv 2 + n 2 q 2 (42) 
with the parameters (cf. the analogous expansion in the 3D case, Ref. (5^|) 



/ 2A(B-E) 
C 2 + 2A(H-F)' 



(43) 



C 2 (B - E) (4A (BH - EF) + C 2 (B + E)) 
V (C 2 + 2A (H - F)f 

The parameter v has the dimensionality of velocity and tends to the first sound velocity in the low-temperature limit 
. As far as the parameter D is proportional to A 2 , the velocity parameter for pair excitations turns to zero at the 
phase boundary when A = 0. In the BEC limit, when Ef, 3> 1, we find that /i — > —Ef,/2 and k —> 1/4. This results 
in the pair excitation spectrum u> q — > q 2 /2 at the BEC side. 

In order to calculate the long-wavelength distribution function, we keep the lowest-order terms in powers of z and 
q 2 in the numerator of Ji w (q, z). This gives us the result 

J to (q,z) = ^^ + -^, (45) 

Z — U! q Z + U! q 

where the coefficients a q and b q are related to the constants determined above as 

_ a + \uj g + xq 2 , _ a - Xuj q + x<? 2 



with the notations 



2lu q [C 2 + 2A (H - F)] ' q 2uj q [C 2 + 2A {H - F)] 



D^D - A^A, \ = A^C~ C^A, X = E^D + D^E - A^B - B^A. (47) 



Here A^, B^, . . . are the derivatives A^ = dA/d/i, etc. The distribution function is calculated setting z = us + id with 
S — > +0. This gives us the structure factor as a superposition of the delta functions. The distribution function then 
takes the form 



4tt C 2 + 2A(H - F) 



« + ^ 2 coth ^ ) _ 1 



Xujq \ 2 



(48) 



For the other distribution function, (q), the derivations are the same, but with a replacement of the derivatives 
over /I by the corresponding derivatives over £. 

C. Distribution functions in the paired state 

Here we consider the paired state of the quasicondensate in which the gap parameter A =/= 0. In this case, the gap 
parameter obeys the gap equation: 

1 r X ( E *l kdk+ l =0 . (49) 



4?r Jo E k g 

The strength g of the contact interaction is expressed through the two-particle binding energy Eb in 2D by the equation 
©. The difference of coefficients A — D is proportional to the left hand side (LHS) of the gap equation. Therefore, 
as long as the gap equation is satisfied, we obtain D = A. Moreover, because the derivatives of matrix elements are 
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calculated while keeping the gap equation satisfied, we find that = and Aq — Z^.This implies in particular 
that the coefficient a — 0, as is evident from expression (|T7)) . and we find: 

Thus gn W ^ (q) tends to a finite value at q — >■ for A ^ 0, and behaves as q~ 2 at q — > for A = 0. As a result, the 
fluctuation contributions nfi and 6n.fi in 2D are finite at A =/= 0. They can diverge only at A = 0. This is to be 
contrasted with the NSR scheme, where the order parameter A is treated as independent variable, and where n/; and 
6n.fi in 2D diverge for all A. 



1 A 

4tt C 2 + 2A (H - F) 




FIG. 1: (Color online) Distribution functions for (a) the fluctuation contribution to the fermion density and (6) the density 
difference, at the binding energy Eb = 0.5-Bf- The solid and dashed curves show the spectral functions obtained, respectively, 
within the GPF formalism and within the standard NSR scheme. The spectral functions are calculated for critical values £ = £ c 
of the chemical potential imbalance, and for two different temperatures. In the graphs, the spectral functions are multiplied by 
q 2 in order to show clearly their behavior at small q. 



Fig. 1 shows the behavior of the fluctuation distribution functions g n (q) and gs n (q) for different temperatures, at 
binding energy Eb — 0.5 and at the critical value of the chemical potential imbalance C = Cc (Eb,T). The critical 
value £ c for a given (Eb, T) is determined as the highest imbalance at which the order parameter A is other than zero. 
The dashed lines correspond to the NSR scheme and reveal a q~ 2 long wavelength divergence. The full lines show the 
results in the GPF scheme, where the long- wavelength divergence is absent. This behavior is seen both for g n (q) and 
gsn(q)- In the limit A — > 0, the functions qg n (q) and qgsn (q) become logarithmically divergent. However, the sign of 
this divergence is opposite to that of the divergence of the functions calculated neglecting the variation of A. 

For the lower temperature shown in Fig. 1, T/Tp = 0.1, g n (q) remains positive, whereas for the higher temperature 
T/Tp = 0.3, there is a sign change in g n (q) as it becomes negative for small q. Regions of negative value for the 
fluctuation distribution function gs n (q) are expected, as the sign will change depending on which species is the majority 
species. However g n (q) is expected to remain positive, as it is proportional to the pair fluctuation density of states. 
The appearance of a long-wavelength instability heralds the breakdown of the paired state. We can track the onset of 
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FIG. 2: (Color online) The fluctuation distribution function g„ (0) as a function of temperature for £ = (a) and for £ = 0.3 
(6). The values of the binding energy are shown in the figure. 



this instability by studying g n (q — > 0) as a function of temperature. At low temperatures, g n (0) is positive, and the 
fluctuation density function remains positive for all q. At high temperatures g n (0) becomes negative, signalling the 
long-wavelength instability. We denote the temperature separating the two regions by T p , and find this temperature 
through solving g n (0) = with respect to temperature. The behavior of <? n (0) as a function of temperature is shown 
in Fig. 2, for different values of the binding energy Ef, and both for balanced and imbalanced systems. The function 
g n (0) diverges when the temperature achieves the limit T — T* at which A = 0. This result explicitly follows 
from the analytic properties of the long-wavelength expansion of the distribution functions as discussed above. The 
temperature T p at which g n (0) = lies below the temperature T* where we find A = 0, and above the critical 
temperature Tbkt for superfluidity. 

The temperature T p does not correspond to a phase transition, because the gap equation is satisfied with a finite 
density both below and above T p . Nevertheless, because T p is the temperature at which the fluctuation density of 
states changes it qualitative behavior, we hypothesize that T p corresponds to a crossover between the normal and 
pseudogap states. This will be further substantiated by comparing our spectral functions to the experimental ones 
in sec. [V] The joint solution of the gap and number equations within the GPF theory then formally provides a 
non-superfluid quasicondensate at temperatures below T p . Indeed, for temperatures Tbkt < T < T p the phase 
coherence is destroyed by the phase fluctuations according to the BKT mechanism, resulting in the phase fluctuating 
quasicondensate discussed by Kagan Q. Through the interpretation of the spectral function, we will denote this 
temperature region as the "pseudogap regime" . It is worth noting that the total fermion density within the GPF 
theory is finite at T = T p without the necessity to introduce any cutoff in the integrals over q. In the next section we 
set up phase diagrams identifying the regions where the superfluid phase and the non-coherent paired phase occur. 
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IV. PHASE DIAGRAMS 



In order to get the complete set of equations for phase diagrams, the number equations ([2"TT) . (|2T|) and the generalized 
gap equation ((5J) are solved jointly with the equation for the BKT transition temperature Tbkt determined by (30l ] 

7T 

Tbkt - -^p s {Tbkt) = 0, (51) 

where p s is the superfluid pair density given by Eq. (|22[) . To investigate the phase transitions for the Fermi gas in 
2D for different binding energies, we have calculated the critical temperatures of the BKT phase transition Tbkt and 
the critical temperature T p below which the phase fluctuating quasicondensate is formed, as a function of the binding 
energy E^. Because the fluctuation contribution to the density is finite at T p and at Tbkt > these temperatures can be 
self-consistently determined from the joint solution of the gap and number equations with the complete thermodynamic 
potential fl = Q sp + Qfi uc t- 
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FIG. 3: (Color online) Phase diagrams for the Fermi gas in 2D (a) in the case of equal spin up and spin down populations, 
(6) for the imbalanced Fermi gas with the chemical potential imbalance £ = 0.5. The crossover pairing temperature T v for 
the pseudogap formation and the BKT transition temperature Tbkt are shown with solid and dot-dashed and solid curves, 
respectively. The dashed and dotted curves show the mean-field phase transition temperatures T* lt T* 2 (explained in the text). 
The arrow indicates the lowest binding energy at nonzero imbalance when preformed pairs can arise. 



The phase diagrams in Fig. 3 show the critical temperatures for cold fermions in 2D as a function of the binding 
energy E\> for the balanced case (panel a) and for the chemical potential imbalance ( = 0.5 (panel b). The formation 
of the superfluid state is indicated by the critical temperature Tbkt of the BKT phase transition. The pseudogap 
temperature T p is the upper bound for the existence of the phase fluctuating quasicondensate described in the previous 
section. We also show the mean- field temperature for pair formation, T*, obtained by solving gap and number 
equations with Vl — Q sp . According to Ref. [54| (for 3D), the unitary gas can exist in the normal state with pairing 
correlations called preformed pairs which survive at temperatures up to this T*. The critical temperatures for the 
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balanced case, Fig. 3 (a), were calculated in Ref. [Hj]. Here, they are reproduced in order to compare them with those 
for a nonzero imbalance. The population imbalance brings new features to the phase diagram: a phase separation 
region (between T*± and T* 2 ) and a minimum binding energy Ef, iCr required for superfluidity. 

For the balanced Fermi gas the superfluid state exists for any value of the binding energy: the BKT critical 
temperature as well as other critical temperatures gradually decrease with decreasing Eb, remaining always finite. 
However, when ( ^ a minimum value of the binding energy Eb,cr is required for superfluidity to exist even at 
T = 0. As shown in Fig. 3, the pseudogap temperature T p does not grow unboundedly when increasing the binding 
energy Ef,. For £ = 0.5 it achieves its maximum at around Eb ~ 4 and then slowly decreases tending to a Unite 
value. Consequently, in the strong-coupling regime the pseudogap temperature is suppressed with respect to the 
mean- field prediction, where it is often identified with our pair formation temperature T*, as in 54j . This behavior is 
qualitatively similar to that for the critical temperature T c as a function of l/a s for the cold fermions in 3D obtained 
first in Ref. Q accounting for the Gaussian fluctuations. 

The critical temperatures T* x and T* 2 coincide with each other in the balanced case, and they can be different 
in the imbalanced case: the area between T* x and T* 2 is the "phase-separated state" . In the phase-separated state, 
uniform phases are not possible. The temperatures T* Y and T* 2 were already calculated in Ref. [33| . The temperature 
T p is determined for the state with A^O. Therefore a non-zero imbalance does not lead to a splitting of this critical 
temperature. However, a tricritical point appears at T p — Tbkt in the phase diagram joining three regions: the 
supernuid state, the pseudogap regime and the normal state. This tricritical point is rather conventional as far as the 
pseudogap temperature indicates a crossover rather than a sharp transition. 




FIG. 4: (Color online) Finite-temperature phase diagram for cold fermions in 2D in the variables (T, Q for the binding energy 
Eb — 0.04i5_F. The full dot indicates a tricritical point. 



At zero imbalance, T p > Tbkt, and the phase coherence in the range Tbkt < T < T p is destroyed by phase 
fluctuations that lead to a phase fluctuating quasicondensate. However, at nonzero imbalance, there is a region where 
pseudogap temperature crosses the BKT temperature for superfluidity. This result is interesting in connection with 
recent experiments on high-T c superconductors (56j , that show a crossing of the zero-field superconducting transition 
temperature and the temperature indicating the opening of the pseudogap in overdoped I^-zSr^CuO^ The crossing 
of pseudogap temperature and BKT temperature is also seen in Fig. 4, showing the phase diagram in the variables 
(T, £), for the binding energy Eb/Ep — 0.04. Here, the same critical temperatures and phase regions are identified as 
in Fig. 3(6). Increasing imbalance is not only detrimental to the superfluid phase, it also suppresses the pseudogap 
regime. 



V. COMPARISON WITH EXPERIMENT 



In the experiment [l7j on pairing of cold fermions in two dimensions, the single-particle spectral function A (q, oj) is 
measured for different values of the wave number q. The spectral function exhibits peaks whose positions indicate the 
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energies of the pair excitations. In the strong-coupling regime, these energies are close to the pair binding energy 
However, as stated in the paper, some discrepancies remain between the peak positions observed in the experiment 
and those predicted by the mean-field theory. The deviation "could stem from beyond mean-field effects provoked by 
our two-dimensional geometry and interaction energy shifts" [l7| . 

In the GPF approach, the pair fluctuation contribution of the fermion density is expressed through the integral 
(|32|) . where the structure factor J (q, uj) describes the spectrum of the pair excitations of the fermion system. Thus 
there should be a correspondence of the peaks of the structure factor J (q, u) with the peaks of the spectral function 
A (q, oj) . In this connection, we compare the positions of the peaks of the spectral function measured in Ref. [l7j 
with those of the structure factor calculated within the GPF approach. The results are shown in Fig. 5 for q = and 
T/T F = 0.27, where k B T F = E F . The 2D scattering length a 2 D is related to the binding energy E b as a 2 D = hj\JmE h . 
The value of the Fermi wave vector taken from Ref. [TtJ is kp = 8.1 |xm" 1 . When using the mass of the fermion atom 
m w 39.964 u, we found that the frequency vp = Epj (27rfi.) corresponding to the Fermi energy is vp = 8.2967 kHz. 

For the visualization of the peaks of the structure factor, we have used J (q, uj + i^) with a finite damping parameter 
7 (as in Refs. [H, |39|, where this parameter was introduced to facilitate the numeric calculations). Here, the value 
7 = 0.2tt/(3 is used, where (3 = 1/ (ksT) is the inverse temperature. 

The parameters of the state (the chemical potential fj, and the gap parameter A) are determined for each plot from 
the joint solution of the gap and number equations. In the number equation, the Gaussian fluctuations are included 
within the GPF formalism. The GPF method provides a finite (convergent) pair fluctuation contribution for any 
finite A without any cutoff for the pair momentum. This is to be contrasted with the standard NSR scheme which 
leads to a divergence of the fluctuation contribution at any A. Therefore the standard NSR scheme cannot be used 
for the description of the pseudogap state, whereas the GPF approach can describe this regime. 
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Frequency (kHz) 



FIG. 5: Full dots: measured energy distribution curve A (q = 0,u) for \n(kFa,2D) 
the fit by elementary functions to the experimental data performed in Ref. [171 ]. 
J (q = 0, ai) calculated in the present work within the GPF approach. 



= 0.8 from Ref. [13]. 
The black solid line: 



The red solid line is 
the structure factor 



In Fig. 5, the high peak at uj = in our results has no relation to the energies of the pair excitations: it is an 
intrinsic feature of the structure factor. The other peak of our structure factor at uj < is positioned remarkably 
close to the measured peak of the spectral function attributed to the pair excitation energy in Ref. [TtJ , especially 
for the relative high coupling strength at In (kpa2D) — 0. A possible reason for the remaining difference between 
the peak positions of the calculated structure factor J (q, uj) and the measured peak spectra can be the experimental 
uncertainty in the determination of the Fermi wave vector, which can be slightly different from the reported value 
kp = 8.1/im _1 . Another possible source of the remaining difference is the similar experimental uncertainty on 
\n{kpa2D)- In particular, this uncertainty can be provided by the facts that the Fermi wave vector determined in 
Ref. [ItJ is a trap-averaged rather than local quantity. It should be noted that the structure factor J (q, w) calculated 
with the mean-field values for /i and A leads to a large discrepancy between the peaks of J (q, ui) and those of the 
measured spectral function. This confirms the importance of including fluctuations through the GPF approach for 
the description of the pseudogap state of cold fermions in 2D. 
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FIG. 6: Solid curve: calculated pseudogap pairing temperature T p (in units of Et/kB) compared with the experimentally [TtJ 
determined pairing temperature T B (full dots). Dotted curve: the mean-field critical temperature T* scaled by the factor 0.36. 



In Ref. [17] . the pairing crossover temperature T* and the pseudogap pairing temperature TJj < T* have been 
introduced. The temperature T* coincides with the mean- field transition temperature T*. The temperature T^, 
as stated in Ref. [17) . indicates the formation of pairs, and has the same physical meaning as the temperature T p 
obtained in our study. As far as the transition between the normal and paired states is a crossover rather than a true 
phase transition, the pairing temperatures Tg and T p only approximately indicate the formation of a paired state. 
In Fig. 6, the pseudogap pairing temperature T p is compared with the experimental data for Tg. The dotted curve 
shows the scaled mean- field transition temperature from Ref. [l7|; the scaling indicates that the experimental result 
for Tg is a factor 0.36 smaller than the mean-field prediction. We see that, in contrast to the mean-field result, the 
value of T p obtained in the present treatment lies in the same range as the experimentally determined temperature 
Tg. This coincidence is worth remarking. However, the conclusions from the latter comparison of two temperatures 
need care, because the temperature in [1 71 ] is measured in the weakly interacting regime and hence it may differ from 
the actual temperature in the strongly interacting regime. 

It is stated in Ref. jlTj that the discrepancy between the mean-field and experimental pairing temperatures could 
suggest that the appearance of a back-bending feature in the spectral function (20| . which has been interpreted 
as a signature for many-body pairing, is only a qualitative evidence. However, the present results show that the 
fluctuations can drastically reduce the pairing temperature T p with respect to T*. Thus there is no discrepancy 
between experiment and theory when taking account the fluctuations. 



VI. CONCLUSIONS 



The T-matrix approach straightforwardly applied to cold fcrmions in two dimensions leads to a divergent fermion 
density for any finite temperature. We have shown in the present work that ta king in to account the variation of 
the order parameter in the number equations, as suggested in the GPF approach [I3l4l5|. provides a divergence- free 
description of the paired state in two dimensions. This was shown both through numerical calculations and through 
an analytic expansion at long wavelengths and low energies, where the divergency occurs in the standard Nozieres 
& Schmitt-Rink approach. The formalism allows to study the effects of the fluctuations both at zero and at finite 
temperatures, and we find that fluctuations affect the critical binding energy to obtain pairing and superfluidity in 
the presence of imbalance. Moreover, the formalism also gives access to the density of states of the pair fluctuations, 
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from which we have defined a pseudogap temperature T p as the temperature where an instability appears in the pair 
fluctuation density. The pseudogap temperature defined in this way agrees with the measured values of the pseudogap 
temperature in 2D fermi gases. Also the location of the peaks in the spectral functions for pair fluctuations is shown 
to agree with the experimental observations. The pseudogap temperature T p , along with the critical temperature 
Tbkt for superfluidity and the pair formation temperature T* , have been calculated as a function of binding energy, 
temperature and imbalance, from which we obtain the phase diagram as shown in Figs. 3 and 4. Whereas in mean- 
field the pseudogap temperature is usually identified with the pair formation temperature, we find that the inclusion 
of fluctuations beyond mean field strongly suppresses the pseudogap temperature with respect to the mean-field pair 
formation temperature. Moreover, in the presence of imbalance, the pseudogap temperature may cross the BKT 
temperature for superfluidity. The results obtained here in the context of supcrfluid quantum gases shed new light on 
the study of the pseudogap phase in layered high-temperature superconductors, where the question of the crossing of 
the pseudogap temperature with the superconducting temperature, and the presence of preformed pairs, remains an 
open question. 
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